##########################
##
##	lab 9 - Simplified SIR model 
##	(should be SEIR model - not all I are 
##	infectious (latent period of 1 week)
##
##	@author: Conrad Stack
##
##########################

##### Problem 1 #####
#Read in measles dataset
measles = read.table("V:/datasets/lab9.dat", sep="\t",header=T)

# line-plot of # cases over time
plot(measles$time, measles$I, xlab="time (years)", ylab="# Infecteds", type="l")

# line-point-plot of # case vs. susceptibles
plot(measles$S, measles$I, xlab="# Susceptibles", ylab="# Infecteds", type="b")

#generate lagged dataset (phase plane)
plot(measles$I[1:545],measles$I[2:546], xlab="I(t-1)", ylab="I(t)", type="b")


##### Problem 2 #####
Inow = measles$I[2:546] 	#I(t)
Sthen = measles$S[1:545] 	#S(t-1)
Ipast = measles$I[1:545]/3300000 #I(t-1)/N

# Inow = beta*Sthen*Ipast = y = (slope)*x

# make a plot - indicates strong positive correlation
SIvect = Sthen*Ipast
plot(SIvect, Inow) 

# find correlation between Inow, SI-vector
corrISI = cor(SIvect, Inow, use="all.obs")
corrISI # 0.97285

# find correlation between Inow, Ipast
# how does it compare to corrISI
corrII = cor(Ipast, Inow, use="all.obs")
corrII # 0.949058

# get regression without intercept
# fit = lm(y~-1+x)
fit = lm(Inow ~ -1 + SIvect)
# slope = transmission rate (B) = 28.1912 ?
B = fit$coef[1]
# R^2 = 0.967

#predict number of cases and plot
# predict(fit)
plot(Inow,predict(fit), xlab="Observed cases", ylab="predicted cases")
# cor.test(Inow,predict(fit), use="all.obs"), should it be exactly the same

##### BONUS #####

# R0 = Inow = BIthenS/N
R0 = 28.19*(3300000-1)*1/3300000
# 28.1899





